# Import and clean incident damage data for a few early San Diego fires
pacman::p_load(sf, tidyverse, fst, stringdist)

rm(list = ls())

source("code/globals.R")

sd_boundary <- st_read(file.path(INPUT_PUBLIC, "CensusTiger/tl_2019_us_county/tl_2019_us_county.shp")) %>%
  filter(GEOID == "06073") %>%
  st_transform(crs = 4326)

sd_damages <- file.path(INPUT_PUBLIC, "DamageData/NotDINS/CAnotDINS/California_SanDiego/FireBurnLayers")


# Paradise (2003) ----
paradise <- st_read(file.path(sd_damages, "Fire_Damage_Paradise.shp")) %>%
  st_transform(crs = 4326) %>%
  dplyr::select(APN, OUT_DM, OUT_DS, SQ_FT, ST_NAME, ST_NUMB, STRUCT_DM, STRUCT_DS, X, Y, REPLACE_CO, CONTENTS_) %>%
  mutate(FireDate = as.Date("2003-10-26"), FireName = "Paradise2003", FIPS = "06073", State = "CA") %>%
  mutate(Address = paste(ST_NUMB, ST_NAME, sep = " ")) %>%
  mutate(NotesField = paste("NumStDest", STRUCT_DS, "ReplaceCost", REPLACE_CO, "ContentsLoss", CONTENTS_, sep = "_"))

# Get geometry and save as lon_orig and lat_orig
paradise <- paradise %>% bind_cols(st_coordinates(paradise) %>% as.data.frame() %>%
  transmute(lon_orig = X, lat_orig = Y))
  
# Make a level of damage variable
paradise$DamageLevel <- NA
paradise$DamageLevel[paradise$STRUCT_DS > 0] <- "Destroyed"
paradise$DamageLevel[paradise$STRUCT_DS == 0 & paradise$STRUCT_DM > 0] <- "Minor"
paradise$DamageLevel <- factor(paradise$DamageLevel, levels = c("Minor", "Destroyed"), ordered = TRUE)

# Limit to parcels with destroyed or damaged main structures, not outbuildings), and trim to final vars of interest
paradise <- filter(paradise, STRUCT_DM > 0 | STRUCT_DS > 0) %>%
  dplyr::select(APN, Address, DamageLevel, FireDate, FireName, FIPS, State, NotesField, lon_orig, lat_orig)


paradise <- group_by(paradise, APN) %>%
  arrange(APN, desc(DamageLevel)) %>%
  mutate(howmanystructures = n()) %>%
  mutate(badrecordflag = howmanystructures > 1) %>%
  filter(row_number() == 1) %>%
  ungroup()

# Witch and associated (2007) -----

witch <- st_read(file.path(sd_damages, "FIRESTORM_2007_DAMAGES_CN.shp")) %>%
  st_transform(crs = 4326) %>%
  mutate(
    SITUS_ADDR = tolower(SITUS_ADDR), SITUS_FRAC = tolower(SITUS_FRAC), SITUS_SUIT = tolower(SITUS_SUIT),
    SITUS_PRE_ = tolower(SITUS_PRE_), SITUS_STRE = tolower(SITUS_STRE), SITUS_SUFF = tolower(SITUS_SUFF)
  ) %>%
  mutate(Address = paste(SITUS_ADDR, SITUS_FRAC, SITUS_SUIT, SITUS_PRE_, SITUS_STRE, SITUS_SUFF, sep = " ")) %>%
  mutate(Address = gsub(" NA ", " ", Address, fixed = TRUE)) %>% # converting everything to lower case before merging makes it easy to pick out the NA values without actually stripping content
  mutate(Address = gsub(" NA ", " ", Address, fixed = TRUE)) %>%
  mutate(Address = gsub(" NA ", " ", Address, fixed = TRUE)) %>%
  mutate(Address = gsub("NA", "", Address, fixed = TRUE)) %>%
  mutate(Address = trimws(Address)) %>%
  rename(APN = APN_8) %>%
  mutate(FireDate = as.Date("2007-10-21"), FireName = "Witch2007", FIPS = "06073", State = "CA") %>%
  mutate(NotesField = COMMENTS)

# Save original lon/lat
witch <- witch %>% bind_cols(st_coordinates(witch) %>% as.data.frame() %>%
  transmute(lon_orig = X, lat_orig = Y))

## - Code the damage amounts
witch$DamageLevel <- NA
witch$DamageLevel[witch$PRI_DES > 0] <- "Destroyed" # See pdf documentation in the data folder - The PRI_DES field lists the number of 'primary structures' destroyed
witch$DamageLevel[witch$PRI_DES == 0 & (witch$PRI_MOD + witch$PRI_MIN > 0)] <- "Minor" # Number of primary structures moderately or minimally damaged
## -- fix some records that appear to have been appended from a different source and have the damage information in the comments field
witch$DamageLevel[witch$COMMENTS == "1 House Destroyed"] <- "Destroyed"
witch$DamageLevel[witch$COMMENTS == "Structure totally destroyed"] <- "Destroyed"
witch$DamageLevel[witch$COMMENTS == "Structure partially destroyed"] <- "Minor"

witch$DamageLevel <- factor(witch$DamageLevel, levels = c("Minor", "Destroyed"), ordered = TRUE)
table(witch$DamageLevel, useNA = "always")

## -- Keep only destroyed and damaged primary homes
witch <- witch[witch$DamageLevel %in% c("Minor", "Destroyed") & !is.na(witch$DamageLevel), ]

## Flag APNs with multiple structures, then keep one record from each group (noting to drop all these later after ZTRAX merge)
witch <- group_by(witch, APN) %>%
  arrange(APN, desc(DamageLevel)) %>%
  mutate(howmanystructures = n()) %>%
  mutate(badrecordflag = howmanystructures > 1) %>%
  filter(row_number() == 1) %>%
  ungroup()

witch <- dplyr::select(witch, APN, Address, DamageLevel, FireDate, FireName, FIPS, State, NotesField, howmanystructures, badrecordflag, lon_orig, lat_orig)

# Add to cleaned dataset and recode some things for consistency with DINS, so that I can process this along with more recent CA fires in subsequent code.
sd <- rbind(paradise, witch) %>%
  rename(UnformattedAssessorParcelNumber = APN) %>% # To match Zillow data naming)
  rename(INCIDENTNAME = FireName) %>%
  rename(date = FireDate) %>%
  mutate(FIPS = as.numeric(FIPS)) %>%
  mutate(INCIDENTNUM = INCIDENTNAME)

rm(paradise, witch)

sd <- sd %>%
  mutate(
    APN_orig = UnformattedAssessorParcelNumber,
    address_orig = Address,
    in_damage_data = 1
  )

sd$DAMAGE <- as.character(sd$DamageLevel)
sd$DAMAGE[sd$DAMAGE == "Destroyed"] <- "Destroyed (>50%)"
sd$DAMAGE[sd$DAMAGE == "Minor"] <- "Minor (10-25%)"
sd$DamageLevel <- NULL

sd$STRUCTURETYPE <- "NotReported"
sd$badstructuretype <- FALSE
sd$ROOFCONSTRUCTION <- "NotReported"
sd$EAVES <- "NotReported"
sd$EXTERIORSIDING <- "NotReported"
sd$WINDOWPANE <- "NotReported"
sd$VENTSCREEN <- "NotReported"
sd$DECKPORCHONGRADE <- "NotReported"
sd$DECKPORCHELEVATED <- "NotReported"
sd$PATIOCOVERCARPORT <- "NotReported"
sd$FENCEATTACHEDTOSTRUCTURE <- "NotReported"

## Save the data
saveRDS(sd, file.path(WORKING, "damage_sd_03_07.RDS"))
